%b00(x)u0''+h00(x)u0'+b10(x)u1''+h10(x)u1'+b20(x)u''+h20u'+b30u''+h30u'=f0
%with Dirichlet boundary 
clear all

%construct mesh
n=400; term=3;
unimesh=(0:1/n:1)';

%construct coefficient function handel  
% 
  thebij;
%  load funcbij;
  thehij;
%  load funchij;
%construct coefficient matrix A
%A=[matelt(b00,h00,n),matelt(b10,h10,n);matelt(b01,h01,n),matelt(b11,h11,n)];
% 
%  A=[matelt(bij{1}{1},hij{1}{1},n), matelt(bij{2}{1},hij{2}{1},n); ...
%      matelt(bij{1}{2},hij{1}{2},n), matelt(bij{2}{2},hij{2}{2},n)];
% 
% A=[ matelt(bij{1}{1},hij{1}{1},n), matelt(bij{2}{1},hij{2}{1},n), matelt(bij{3}{1},hij{3}{1},n), matelt(bij{4}{1},hij{4}{1},n); ...
%     matelt(bij{1}{2},hij{1}{2},n), matelt(bij{2}{2},hij{2}{2},n), matelt(bij{3}{2},hij{3}{2},n), matelt(bij{4}{2},hij{4}{2},n);...
%     matelt(bij{1}{3},hij{1}{3},n), matelt(bij{2}{3},hij{2}{3},n), matelt(bij{3}{3},hij{3}{3},n), matelt(bij{4}{3},hij{4}{3},n);...
%     matelt(bij{1}{4},hij{1}{4},n), matelt(bij{2}{4},hij{2}{4},n), matelt(bij{3}{4},hij{3}{4},n), matelt(bij{4}{4},hij{4}{4},n)];

A=matA(n,term);
if(size(A,1)==0)
    disp('too high')
else
    

%construct source f
b=srcf(n,term);

%solve the equation
u=A\b;

%measure the error
if(term==4)
uex=[u0ex(unimesh(2:n));u1ex(unimesh(2:n));u2ex(unimesh(2:n));u3ex(unimesh(2:n))];
elseif(term==1)
    uex=u0ex(unimesh(2:n));
elseif(term==2)
    uex=[u0ex(unimesh(2:n));u1ex(unimesh(2:n))];
elseif(term==3)
    uex=[u0ex(unimesh(2:n));u1ex(unimesh(2:n));u2ex(unimesh(2:n))];
end
err=max(abs(uex-u))


%uex=uexact(n,term);
%err=L2error(uex,u)

%  err=maxL2err(u,term)
end